/*------------------------------------------------------------------------
 *  Copyright (C) 2011  Luis M. de la Cruz Salas
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

#ifndef _MOMENTUM_XCODI_H_
#define _MOMENTUM_XCODI_H_

#include "Equations/GeneralEquation.hpp"

namespace Tuna {

template <typename T>  class Momentum_XCoDi;
  
template <typename T>  
struct Typeinfo <Momentum_XCoDi<T> > {
  typedef typename Typeinfo<T>::prec_t prec_t;
  enum { Dim = Typeinfo<T>::Dim };
};

/*! 
 ***************************************************************************
 *  This class derives from GeneralEquation through B & N trick.
 *  In this class some particular variables and methods for scalar PDE's
 *  are defined. This class implements the FVM discretized x-component of 
 *  momentum equation for Newtonian incompressible fluids.
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Thu May  7 13:55:49 CDT 2009 ]
 ***************************************************************************
 */
template<typename Tscheme>
class Momentum_XCoDi : public GeneralEquation<Momentum_XCoDi<Tscheme> >
{
  using GeneralEquation<Momentum_XCoDi<Tscheme> >::name;
    
public:
  typedef typename Typeinfo<Tscheme>::prec_t prec_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::tiny floatTinyArray_t; 
  typedef typename TunaArray<int, Typeinfo<Tscheme>::Dim >::tiny intTinyArray_t;
  typedef typename TunaArray<prec_t, Typeinfo<Tscheme>::Dim >::huge ScalarField;
  typedef SparseMatrix< Diagonal< prec_t, Typeinfo<Tscheme>::Dim > > DiagMat;

  Momentum_XCoDi() : GeneralEquation<Momentum_XCoDi<Tscheme> >() { }

/*!
 *  Main constructor.
 *  The scalarfield \c phi_0 defined here is a reference to a global field 
 *  variable. Here \c phi_0 is used as an old value of that field variable. 
 */
  Momentum_XCoDi(ScalarField &phi_global,
		       DiagMat &matrix, ScalarField &b,
		       const floatTinyArray_t& deltas)
    : GeneralEquation<Momentum_XCoDi<Tscheme> >(phi_global,matrix,b,deltas)
  {
    Dim = Typeinfo<Tscheme>::Dim;
    phi_0.reference(phi_global);
    name = "Momentum_XCoDi";
    Gamma = 1.0;
  }

  ~Momentum_XCoDi() { }; // Nothing to do  
/*!
 *  The Curiously Recursive Template Pattern (CRTP) is used.  
 *  asDerived() gets an returns a reference to the derived class.
 */
  inline Tscheme& asDerived() { return static_cast<Tscheme&>(*this); }
  
//  
// Set functions used to initialize the parameters.
//
  inline void setGamma(prec_t G) {Gamma = G;}
  inline void setRayleigh(prec_t R){Rayleigh = R;}
  inline void setUvelocity(const ScalarField& uu) { u.reference(uu); }
  inline void setVvelocity(const ScalarField& vv) { v.reference(vv); }
  inline void setWvelocity(const ScalarField& ww) { w.reference(ww); }
  inline void setTemperature(const ScalarField& TT) { T.reference(TT); }
  inline void setPressure(const ScalarField& pp) {
    du.resize(pp.shape());
    p.reference(pp);
  }  
//
// Get functions
//
  inline ScalarField& get_du() { return du; }
/*!
 *  Call to the functions defined in the derived classes.
 */
  inline bool calcCoefficients() 
  {
    if(Dim == 1)
      return asDerived().calcCoefficients1D();
    else if (Dim == 2)
      return asDerived().calcCoefficients2D();
    else if (Dim == 3)
      return asDerived().calcCoefficients3D();
  }

///
/// Function to calculate the coefficients of discrete equations for
/// turbulent regime. Turbulent viscosity is passed as argument (nut)
///
  inline bool calcCoefficients(const ScalarField &nut) {
    return asDerived().calcCoefficients(nut);
  }   

  inline void calc_du_2D( );
  inline void calc_du_3D( );
  /****
  void NUcalc_du_2D( );
  void NUcalc_du_3D( ); 
  ****/

  inline void printInfo() {
    std::cout << "\n +-----+"
	      << "\n | Momentum_XCoDi< ";
    asDerived().printInfo();
    std::cout << "< " << Dim << "D  > >"
	      << "\n +-----+"
	      << "\n | Gamma    = " << Gamma
	      << "\n | Rayleigh = " << Rayleigh
	      << "\n +-----+";  
  }

protected:    
  int Dim;
  prec_t Gamma, Rayleigh;
  ScalarField T, u, v, w, p, du;
  ScalarField phi_0;

};

template<typename Tscheme>
void Momentum_XCoDi<Tscheme>::calc_du_2D( ) {
  Range I(this->bi, this->ei), J(this->bj, this->ej);
  
  du(I,J) = this->dy / ( this->aP(I,J) - (this->aE(I,J) +
					  this->aW(I,J) +
					  this->aN(I,J) +
					  this->aS(I,J)) );
}

template<typename Tscheme>
void Momentum_XCoDi<Tscheme>::calc_du_3D( ) {
  prec_t dyz = (this->dy) * (this->dz);
  Range I(this->bi, this->ei), J(this->bj, this->ej), K(this->bk, this->ek);
    
    du(I,J,K) = dyz / ( this->aP(I,J,K) - (this->aE(I,J,K) +
					   this->aW(I,J,K) +
					   this->aN(I,J,K) +
					   this->aS(I,J,K) +
					   this->aF(I,J,K) +
					   this->aB(I,J,K)) );
}
  
/**** THIS WILL BE REVISED TO MANAGE NON-UNIFORM MESHES
template<typename Tscheme>
void Momentum_XCoDi<Tscheme>::NUcalc_du_2D( ) {
//	Range I(this->bi, this->ei),
//	    J(this->bj, this->ej);
    for(int j = this->bj; j <= this->ej; ++j)
	for(int i = this->bi; i <= this->ei; ++i) {
	    du(i,j) = this->dyn(j) / ( this->aP(i,j) - (this->aE(i,j) +
							this->aW(i,j) +
							this->aN(i,j) +
							this->aS(i,j)) );
	}
}

template<typename Tscheme>
void Momentum_XCoDi<Tscheme>::NUcalc_du_3D( ) {
//	Range I(this->bi, this->ei),
//	    J(this->bj, this->ej),
//	    K(this->bk, this->ek);
    for (int k = this->bk; k <= this->ek; ++k)
	for (int i =  this->bi; i <= this->ei; ++i)
	    for (int j = this->bj; j <= this->ej; ++j) {
		du(i,j,k) = (this->dyn(j)) * (this->dzn(k)) /
		    ( this->aP(i,j,k) - (this->aE(i,j,k) +
					 this->aW(i,j,k) +
					 this->aN(i,j,k) +
					 this->aS(i,j,k) +
					 this->aF(i,j,k) +
					 this->aB(i,j,k)) );
	    }
}
****/

}  // namespace Tuna


#endif //_MOMENTUM_XCODI_H_













